#!/usr/bin/perl
# 2012/11/12 removing potentially PCR biassed data (if linker is different, allow duplicate)

use strict;
use warnings;
use IO::File;
use Getopt::Std;
use Carp qw(croak);
$| = 0;

if(@ARGV != 4 or $ARGV[0] eq '--help'){
	die "Usage : $0 -i [map file] -o [output filtered map file]\n";
}

my %opt;
getopts("i:o:", \%opt);
my $FILE_map = $opt{i};
my $FILE_out = $opt{o};


my $COUNTER_BIAS = 0;
my $COUNTER_UNIQUE = 0;

my $before = '';
my %foundLinker;
my $fh_out = IO::File->new($FILE_out, 'w') or die "cannot write $FILE_out: $!";
my $fh_map = IO::File->new($FILE_map) or die "cannot open $FILE_map: $!";
while($_ = $fh_map->getline()){
#	s/\r?\n//;	# 改行を加えるためにわざとコメントアウト
	my ($id, $chr1, $loc1, $direction1, $mapQ1, $hinID1, $hinLoc1, $chr2, $loc2, $direction2, $mapQ2, $hinID2, $hinLoc2, $type, $num, $locations) = split /\t/;

	# both side should have data
	if($chr1 eq 'NA' or $chr2 eq 'NA'){
		next;
	}

	# check id
	my $linker = "";
	if($id =~ m/ChIA-PET_linker(\w\w)/){
		$linker = $1;
	}


	my $current = join("\t", $chr1, $loc1, $chr2, $loc2);
	if($before ne $current){
		$fh_out->print("$_");
		%foundLinker = ();
		$COUNTER_UNIQUE++;
	}else{
		if(exists $foundLinker{$linker}){
			$COUNTER_BIAS++;
		}else{
			$fh_out->print("$_");
			$COUNTER_UNIQUE++;
		}
	}
	$before = $current;
	$foundLinker{$linker} = 1;
}

$fh_out->close();
$fh_map->close();

my $SUM = $COUNTER_UNIQUE + $COUNTER_BIAS;
print "PCR biasness remove\n";
printf "Unique: %d (%.3f)\n", $COUNTER_UNIQUE, $COUNTER_UNIQUE / $SUM * 100;
printf "Bias: %d (%.3f)\n\n", $COUNTER_BIAS, $COUNTER_BIAS / $SUM * 100;

